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ABSTRACT 

We describe the screened Korringa-Kohn-Rostoker (KKR) method and the third- 
generation linear muffin-tin orbital (LMTO) method for solving the single-particle 
Schrodinger equation for a MT potential. In the screened KKR method, the eigenvec- 
tors crl i are given as the non-zero solutions, and the energies £j as those for which such 
solutions can be found, of the linear homogeneous equations: J2rl Kr>l',rl ( £ i) c RL,i = 0, 
where K a (e) is the screened KKR matrix. The screening is specified by the boundary con- 
dition that, when a screened spherical wave ip RL (e, r R ) is expanded in spherical harmonics 
Yr'u Otr') about its neighboring sites R', then each component either vanishes at a radius, 
r R l=( iR'L'i or is a regular solution at that site. When the corresponding "hard" spheres are 
chosen to be nearly touching, then the KKR matrix is usually short ranged and its energy 
^sO \ dependence smooth over a range of order 1 Ry around the centre of the valence band. The 
KKR matrix, K (e u ) , at a fixed, arbitrary energy turns out to be the negative of the Hamil- 
tonian, and its first energy derivative, K (e u ) , to be the overlap matrix in a basis of kinked 
partial waves, (zv, i\r) , each of which is a partial wave inside the MT-sphere, tailed with 
a screened spherical wave in the interstitial, or taking the other point of view, a screened 
spherical wave in the interstitial, augmented by a partial wave inside the sphere. When of 
short range, K (e) has the two-centre tight-binding (TB) form and can be generated in real 
space, simply by inversion of a positive definite matrix for a cluster. The LMTOs, xrl [sv) > 
are smooth orbitals constructed from § RL (e„, r^) and $^ (e„, r R ) , and the Hamiltonian 
and overlap matrices in the basis of LMTOs are expressed solely in terms of K (e„) and its 
first three energy derivatives. The errors of the single-particle energies £, obtained from the 
Hamiltonian and overlap matrices in the $ (e v )- and x { £ v) bases are respectively of second 
and fourth order in £j — e v . Third-generation LMTO sets give wave functions which are 
correct to order — e u , not only inside the MT spheres, but also in the interstitial region. 
As a consequence, the simple and popular formalism which previously resulted from the 
atomic-spheres approximation (ASA) now holds in general, that is, it includes downfolding 
and the combined correction. Downfolding to few-orbital, possibly short-ranged, low-energy, 
and possibly orthonormal Hamiltonians now works exceedingly well, as is demonstrated for 
a high-temperature superconductor. First-principles sp 3 and sp 3 d 5 TB Hamiltonians for the 
valence and lowest conduction bands of silicon are derived. Finally, we prove that the new 
method treats overlap of the potential wells correctly to leading order and we demonstrate 
how this can be exploited to get rid of the empty spheres in the diamond structure. 
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INTRODUCTION 

There is a need for an intelligible and accurate first-principles electronic-structure method. 
Our efforts have been directed towards developing a single-particle basis which, for the 
application at hand, can be adjusted to a useful compromise between being short ranged, 
minimal, and accurate. 

Recent developments of Multiple Scattering Theory 

Since atoms are nearly round, it seems most natural to start out using spherical waves, 
ji (kx) Yl (?) and n\ (nr) Yl (?) , where L—lm, as done when solving Schrodinger's equation 
with the classical multiple-scattering method due to Korringa, Kohn, and Rostoker (KKR).0 
In this method the scattering by the atom at site R is specified by the phase shifts, t]ri (k) , 
of its potential well, and the structure of the solid is specified by a Hermitian matrix with ele- 
ments Brl^li (/c) , the structure constants. In terms of these, the wave-function coefficients, 
crl,z, are the solutions of the homogeneous, linear equations, one for each R'L' : 

i B R'L',RL («) + « COt T] m («) 5r'L',Rl] C RLii = 0, (1) 

RL 

and the energies, £j, are the values of k 2 (= e) for which solutions can be found, i.e. the 
determinant of B (k) + k cot rj («) vanishes. There are merely 4-16 equations per atom 
because all phase shifts with I > 1-3 vanish. The KKR equations provide the exact solutions 
of Schrodinger's equation, but only for a muffin-tin (MT) potential, V (r) = J2rVr (|r — R|) , 
which is a superposition of spherically symmetric, non-overlapping potential wells, Vr (r) , 
of ranges sr. 

There are three problems with the KKR method: First of all, a non-overlapping MT 
potential is a poor representation of the self-consistent potential in any, except the most close 
packed solid. Secondly, the structure constants have long range, and thirdly, strong energy 
dependence. Specifically, the energy dependence of the KKR matrix, B (k) + K cot 77 (ft) , 
is not linear like that of the secular matrix, H — eO, for an energy eigenvalue problem in 
an energy- independent, possibly non-orthonormal representation. The main reason for the 
second and third drawbacks is that the spherical Bessel (ji) and Neumann (jii) functions 
have long range and depend on energy. This leads to interferences, which cause long range 
and strong energy dependence of the structure matrix. For a crystal with lattice translations 
T, the Bloch-summed structure matrix, Br/l>,rl (k, k) = Z^T ex P(^k-T) -B(_r'+t)l',rl ( k ) , 
must be evaluated by the Ewald procedure, has poles at the free-electron parabola, n 2 = 
J2g |k + G| , and a branch cut at the bottom of the continuum, k—0. 

It was recently shown,! and we shall present a slightly different proof below, that even 
when the potential wells overlap, the KKR equations do hold to first order in the potential 
overlap. This, as we shall demonstrate, allows the use of MT spheres with up to at least 50 
per cent radial overlap [sr + Sri < 1.5 |R — R'| for R and R' denoting nearest neighbors], 
and hence treat the potential between the atoms in a more realistic way. With such large 
overlaps, the zero of the potential moves from the potential threshold between the atoms 
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towards the vacuum level, and this means that the energies for the occupied states are usually 
negative. 

It was furthermore shown,i that transformation to linear combinations of spherical Bessel 
and Neumann functions, so-called screened spherical waves, characterized by a set of back- 
ground phase shifts am (k) , can remove the long range and the strong energy dependence 
from the structure matrix, provided that the energy is not too high. This screening trans- 
formation may be expressed as: 

tana „„\ , . , r „„,_i tana 



\n a ) = \n) 1 B a , where tan 77° = tanr?-tana, and [B a p = B' 1 + . (2) 

V K J K 

These are, respectively, vector-, scalar-, and matrix equations. The superscript labels the 
representation. In the first equation, we have used a notation in which \n) is a row vector 
of functions with components n\ (kxr) Y im (r R ) and where i\r = r — R. The last equation 
involves inversion of the matrix B + Kcota. For a set of background phase shifts, which are 
known to give short range, this inversion may be performed in real space and the screened 
KKR method is basically a first principles tight-binding (TB) method. At the time,i how- 
ever, the relation between range and background phase shifts was poorly understood. This 
problem was solved later! by expressing the background phase shifts in terms of their hard- 
sphere radii, a R i, defined by 

tan a RL («) = ji {na RL ) / m {Ka RL ) , (3) 

or equivalently, by letting the background phase shifts be those of repulsive potential wells.! 
Looked upon in this way, the role of the confinement is to push the bottom of the continuum 
up in energy with respect to the floor of the MT potential, and thereby, to leave below a 
range of energies in which the confined wave-equation solutions are localized and which, in 
order to be useful, should include the range of the occupied bands for the real (attractive) 
potential. With the definition (0), the screened spherical waves, \n a ) , are still quite energy 
dependent, but only due to their normalization. A more suitable normalization followed 
naturally from the hard-sphere point of view.i 

Finally, it was pointed out! that the screening transformation may also be used to remove 
unwanted channels from the KKR equations by choosing, for those channels, the background 
phase shifts equal to the real phase shifts, a (k) = r] (/«). This is a transformation to a minimal 
basis. 

With these three recent developments, we have the basic ingredients for the Schrodinger 
part of an intelligible and accurate method, the third-generation LMTO method. 

Earlier developments; the Atomic Spheres Approximation 

The attempts to develop from the KKR method an intelligible first-principles method 
were initiated 25 years ago! and overlapping spheres, the two-centre interpretation, and 
screening transformations! have been used routinely for a long time. The new development,! 
which started six years ago,! and which will be further elaborated on in the present paper, 
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aims at making the method also accurate without loss of simplicity and elegance. Whereas 
in the earlier developments reduced range and energy dependence were achieved through a 
physically motivated approximation, namely the atomic-spheres approximation (ASA), this 
is now achieved exactly, and that in turn, allows a controlled approximation for the potential 
overlap. 

The ASAii consists of letting the MT spheres overlap to the extent that they become 
space filling, whereby the interstitial region is effectively eliminated so that one may neglect 
the energy dependence of the wave functions in this region and, hence, the energy depen- 
dence of the structure matrix. The remaining energy dependence now occurs only along the 
diagonal of the KKR-ASA matrix where it enters through the radial logarithmic derivatives 
evaluated at the atomic sphere, 

D {(f) Rl (e, s R ) } = s R <p' Rl (e, s R ) / <p Rl (e, s R ) . (4) 

If, as is usually done, the kinetic energy in the interstitial is taken to be zero (k 2 = 0), the 
suitably renormalized KKR-ASA equations become: 

53 [ S R'l>,rl ~ P Ri (e) 5r'l>,rl] c RL ,i = 0, where 



RL 



P% (e) = 2 (21 + 1) (-) 

\Sr/ 



wy^D&Ri ( £ , SR )} + 1 + 1 



+ 1RI 



e-C 



Rl 



,s R J D {(f) m (e,s R )} - I 
are the potential functions for well v R (r) , and C R i, A ra , and 7^ are potential parameters. 
S° is the structure matrix given by: 

S° nsi te = 0, S° SCT = -2 (w/d) , S° spa = 2V3 (w/df , , S° dd{a ^ 5) = 10 (w/df (-6, 4, -1) , (5) 

when we choose R' — R = id and w is an arbitrary length scale, usually chosen to be 
the average Wigner-Seitz radius. The structure matrix thus consists of effective hopping 
integrals. For monatomic crystals, this gave rise to the concept of canonical bands.! However, 
the d^ 1 ^ 1 '^ 1 -decay of the hopping integral between orbitals of angular-momentum characters 
I and I' is too slow for a tight-binding scheme, except for d- and /-orbitals. 

It was therefore a breakthrough when it became understood that similarity transfor- 
mations could be performed on the KKR-ASA (k 2 = 0) equations and could lead to short 
range.! The transformation from the bare to a screened representation, specified by screening 
constants a R L, is given by: 

P a (£)- X = P°(e) -1 -a and [ST 1 = [S°] ~' - a, (6) 

which are respectively scalar- and matrix equations, a is a diagonal matrix with elements 
ctRi. The corresponding transformation for the resolvent, useful for Green-function and CPA 
calculations,! is: 



P b (z) 



PTO PTO rpa ( , _ ST 1 PMfO (7) 
3) Pb (z) + P» (z) [M {Z) bJ Pb(*)' {1) 
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This involves no matrix multiplications, but merely rescaling of matrix elements. As for 
TB theory, the screening transformation gave a formalism for the "kinetic" part of the often 
observed dependence of the hopping integrals and on-site elements on the environment J The 
screening constants yielding short range were found empirically The potential-dependent 
choice ara=7ra, on the other hand, makes the energy dependence of the KKR-ASA ma- 
trix linear (to second order) so that C + a/AS 7 v / A = h 1 becomes the Hamiltonian in an 
orthonormal, but not necessarily short-ranged basis. The transformation finally made it pos- 
sible to remove channels from the KKR-ASA equations by choosing (e) = P° Rl [e)~ for 
such channels.^ This removal, or " downfolding" , however, reintroduced energy dependence 
of the structure matrix. 

When performing density-functional calculations one needs to solve not only 
Schrodinger's but also Poisson's equation, and with the ASA method this involves approx- 
imating not only the potential but also the charge density by a superposition of slightly 
overlapping, spherically symmetric contributions. This gives a very simple scheme which 
fails badly in describing total-energy changes caused by symmetry-lowering distortions,^ 
however, e.g. the ASA can be used for calculation of pressure-volume relationsjl! but not 
for calculation of phonon frequencies. Moreover, since the potential spheres are supposed to 
be space filling in the ASA, open structures can only be treated if the interstices between 
the atoms are filled with "empty" spheres and this works well only for structures such as the 
diamond structure, where the interstices have high symmetry. Even in such a case, for the 
description to be intelligible all empty-sphere channels must be downfolded, and that intro- 
duces a rather strong, non-linear energy dependence of the structure matrix which cannot 
be treated in the LMTO-ASA approach to be discussed below.0 

Linear Muffin-Tin Orbitals of the first and second generations 

In practice one does not solve the KKR equations, but one uses Green functions in 
a short-ranged representation and at complex energies,! or one solves energy eigenvalue 
equations, J2rl[Hr>l',rl — £Or'l',rl}crl,i = 0, which are equivalent with the KKR equa- 
tions in a certain energy range around some chosen energy, e v . In the linear muffin-tin 
orbital (LMTO) methodJlH such an eigenvalue problem is arrived at by using the Raleigh- 
Ritz variational principle for the Hamiltonian in a basis of LMTO's constructed from the 
radial Schrodinger-equation solutions, (f>jy(e,r), for the potential wells and their first en- 
ergy derivatives, 4>m (s,r) , at the chosen energy, e—e v . In the interstitial region, the first 
and second generation LMTOs use the spherical waves at but not their first energy 
derivatives, so that the energy dependence in the interstitial is suppressed. The second- 
generation LMTO formalism! is elegant, but only in the ASA and only if no channels have 
been downfolded. Under these conditions, the Hamiltonian and overlap matrices are ex- 
pressed solely in terms of the structure matrix and the potential functions: The struc- 
ture matrix enters the formalism in the form of a first-order, two-centre TB-Hamiltonian: 
h = P -1 / 2 (S — P) P -1 / 2 , where as usual P (e) is a diagonal matrix. Here and in the follow- 
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ing, the common superscript a is dropped and an omitted energy argument means that the 
energy is set to e v . In terms of this two-centre Hamiltonian the LMTO set may be expressed 
as \x) = \4>) + 0) h, a form which may be regarded as the matrix equivalent of the lin- 
ear approximation (f>m (e, r) rs 4>m ( r ) + 0ra ( r ) ( £ — £ v) ■ In the basis of these LMTOs the 
Hamiltonian and overlap matrices are respectively: 

(x|-A + V-e v \x) = h(l + oh) and (x \ x) = (1 + ho) (1 + oh) + hph, (8) 

where 

o=(0|^)=lp/P and p+o 2 E(^=ip/P (9) 

are diagonal matrices and it has been assumed that (f>m (r R ) Yi m (i R ) is normalized to unity in 
its sphere, i.e. that (0 2 ) = 1. The overlap matrix is seen to be nearly factorized and one may 
therefore transform to the Lowdin-orthonormalized representation, ^ _L \ = \x) (x I x)~^ 2 ~ 
\x) (1 + oh)^ 1 = |x 7 ) , in which one finds the following expansion for the Hamiltonian: 

X L I -A + V - e v \ x L ) = h - hoh + h [oho - (ph + hp) /2] h + ... . (10) 

When a gives short range this is a power series in a TB Hamiltonian, h 3 . Truncation of this 
series after the first term yields a spectrum which is accurate in an energy window of size 
~ (IO0)" 1 = |P/P around e v , but distorted further away. Adding terms, increases the size 
of this window at the expense of including further hoppings. The form (JlOD has been useful 
in recursion calculations!!! for structurally disordered condensed matter.ll3 

For a case like the diamond structure, where one only wants LMTOs centered on atoms, 
downfolding of the empty-sphere LMTOs is achieved by transformation of the structure ma- 
trix using: 3e = P% i e v)~ l > with E referring to the empty-sphere channels. The energy is here 
set to e v because in the LMTO-ASA formalism the structure matrix must be energy indepen- 
dent. Now, an atom-centered LMTO has a tail which extends into the empty spheres, and 
here, it is substituted by the corresponding partial waves. The atom-centered LMTO is there- 
fore: \xa) = \4>a) + <Pa) h AA + \4> E ) h EA , with h EA = [-dP° E /de\ Ev 1/2 S EA (e„) P A 1/2 . 
This is the way in which the energy dependence of a E (e) enters, but only to linear order. The 
overlap matrix {xa I Xa) will now contain the term h AE h EA involving A — E — A hoppings, 
in addition to the terms in (^). This is clumsy and ruins the near factorization of the overlap 
matrix. With downfolding, the power-series expression fllpp for the LMTO Hamiltonian in 
the Lowdin-orthonormalized basis does therefore not apply. 

Most LMTO calculations include non-ASA corrections to the Hamiltonian and overlap 
matrices, such as the combined correction for the neglected integrals over the interstitial 
region and the neglected partial waves of high /. This brings in the first energy derivative of 
the structure matrix, S, in a way which makes the formalism clumsy.!'! Our current, second- 
ereneration LMTO co de§ is useful and quite accurate for calculating energy bands because it 
includes downfolding in addition to the combined correction,0 but the underlying formalism 
is so complicated that we never tried to publish it. On the other hand, the combined 
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correction is often important, and so is downfolding because it is the only accurate means of 
avoiding "ghost bands". The reason for the lost elegance beyond the ASA is that, whereas 
the LMTO basis is complete to first order in e — e v inside the spheres, it is only complete 
to zeroth order in the interstitial. A compact formalism is therefore obtained only when 
the interstitial region is neglected, and that is what the ASA does, simply by substituting 
the MT spheres by space-filling spheres and neglecting the overlap errors. The proof that 
the KKR equations hold to leading order for overlapping potentials^ does not apply to the 
LMTO-ASA formalism. 

There are LMTO methods sufficiently accurate to provide ab initio structural energies 
and forces within density-functional theory.0 For the reason mentioned above, the LMTOs 
for such methods^ are defined with respect to non-overlapping potentials, and since there 
is considerable probability that a valence or conduction electron is in the interstitial region, 
outside atom-centered, non-overlapping spheres, an accurate basis has to include extra de- 
grees of freedom to describe this region, empty-sphere orbitals centered at interstitial sites 
and/or atom-centered LMTOs with tails of different kinetic energies (multiple kappa -sets). 
Moreover, these methods do not use small and short-ranged representations. Finally, since a 
non-overlapping MT potential is a poor approximation to the self-consistent potential, these 
methods must include the matrix elements of the full potential. Hence, the formalisms are 
set up to provide final, numerical results and by themselves provide little insight. 

Third-generation LMTOs 

In this paper we shall modify the LMTO set without increasing its size, in such a way that 
it becomes complete to first order in the interstitial region too. This is a rather natural thing 
to do, once the screened spherical waves have been defined in terms of hard-sphere radii. For 
the MT Hamiltonian, including downfolding, we shall regain the simple formulas from the 
ASA, provided that \<p) , h, o, and p are suitably redefined. The Hamiltonian and overlap 
matrices are now given solely in terms of the screened and renormalized KKR matrix, which 
we shall name K (e„) , and its first three energy derivatives, K (e v ) , K (e„) , and K (s v ) ; the 
potential parameters and the structure matrix do not occur individually as in the formalisms 
of the previous generations. Third-generation LMTOsl thus do satisfy the definition that 
they form a basis constructed to reproduce the wave functions, ^ (r) , for a MT potential 
to linear order in the deviation of the single-particle energy, £j, from a freely chosen level, 
e v . That is, the error of the wave functions is of second order in e$ — e v and the error of 
the single-particle energy is then of fourth order. When we use potential wells that overlap, 
the wave functions will be correct to linear order in the potential overlap and the energy 
error will be of second order. As we shall demonstrate, this will remedy all shortcomings 
mentioned above for the previous LMTO generations. 

We shall only be concerned with solving Schrodinger's equation in the present paper and 
leave our LMTO-like expansion of the charge density, solution of Poisson's equation, and 
evaluation of the total energy and forces for future papersB 
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We start with a concise yet self-contained derivation of the screened KKR method, which 
will lead to suitably renormalized versions of Eq.s (0). Then we derive an expression for the 
error caused by using this method for potential wells which overlap, and find that the error is 
of second order in the overlap. The weak energy dependence and short range of the screened 
and renormalized KKR matrix, K (e) , is exploited by using it to generate few-orbital, low- 
energy, possibly orthonormal and short-ranged Hamiltonians for a generic high-temperature 
superconductor (HTSC). Thereafter we derive the new LMTO method and demonstrate by 
application to free electrons that its energy errors are really of fourth order in — e v . The 
power and flexibility of the new method is demonstrated by deriving for the HTSC and for 
diamond-structured silicon various LMTO sets. Using non-orthogonal sp 3 sets for Si, we can 
get an accurate first-principles description of the valence and conduction bands if a 12th- 
nearest-neighbor range is allowed in the Hamiltonian and overlap matrices. With e v chosen 
in the middle of the valence band, a 6th-n.n. sp 3 -set suffices for an accurate description of 
the valence band and a reasonable description of the conduction band. In order to halve the 
number of matrix elements, even for a non-orthogonal basis, we use a formalism analogous 
to (H) where the off-diagonal elements of o and p have been neglected so that h is the only 
matrix. With this simplification of the Hamiltonian and overlap matrices, retaining the 
6th-n.n. sp 3 basis and the low e u , the description of the valence band remains good and 
merely the conduction band deteriorates. Finally, it is possible to limit the range to 3rd 
nearest neighbors provided that (i-orbitals are included in the basis. In the last section, we 
demonstrate that not only for the KKR method, but also for the new LMTO method, the 
overlap error is of second order and that this can be exploited to get completely rid of the 
empty-sphere wells in the diamond structure. 

SCREENED SPHERICAL WAVES 

We start by defining sets of solutions of the wave equation, [A + e] ij) (e, r) = 0, so-called 
screened-spherical-wave sets, {ipR L (e, r — R)} , which will serve as interstitial (envelope) 
functions for the basis that we shall use for solving Schrodinger's equation. The members 
of a screened-spherical-wave (SSW) set are obtained by letting R run over all atomic sites 
and L over all angular-momenta for which the scattering is strong. The set is labelled by 
the superscript a. Instead of defining %p RL (e, r#) as a specific linear combination of spherical 
Neumann functions like in Eq. (0), we specify it in terms of an inhomogeneous boundary 
condition which is illustrated in Fig.s |l] and |2| and is given as follows: 

Concentric with each MT sphere, R', we imagine a series of possibly coinciding "hard" 
spheres with radii Or/^/. Now, ip RL (e, tr) is that solution of the wave equation whose 
Yr'l 1 projection on the R'L' sphere equals Srl,r>l', that is, 1 on its own sphere and 
on all other spheres. We do not associate SSWs and hard spheres with weakly- and non- 
scattering channels. For such a channel, the Yr/z/ (?r/) projection of the SSWs is defined to 
be a regular solution of the corresponding radial Schrodinger equation, that is, it matches 
onto the irregular wave-equation solution jv (wb') — tan^R^/ (k) riy (ktri) , times some con- 
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stant, c R , L , RL (e) . The weakly- and non-scattering channels are thus parts of the SSW and 
will not enter the screened KKR- and LMTO matrices explicitly. All high-/' channels are non- 
scatterers [tan t)r>i' (k) = 0] due to the dominance of the centrifugal barrier. Empty spheres 
are examples of a weak scatterers. Strong scatterers are then, by definition, those channels 
with which we associate SSWs and hard spheres. Note that all SSWs in the set have the 
same boundary condition, except for the 5rl,r>v- The SSW set, {ip RL (s, tr)} = \RL) (r| , 
may thus be considered as an unperturbed Green function in a hybrid representation. 

Fig. ^ shows an SSW for the hypothetical case of only strong scattering. Weak- and non- 
scattering channels would have shown up as little tails extending into the two hard spheres. 
Such tails may be seen in Fig. ^ where the dashed curve is an SSW for Si. In this figure 
we have set the radial functions of the strongly scattering channels to zero inside the hard 
spheres and, defined in this way, ip RL (e, r R ) jumps by the amount Y L (r R ) at its own hard 
sphere, r R =a R ^, and has kinks at all hard spheres. Had we instead chosen to continue also 
the strongly scattering channels of the SSW into the hard spheres, the SSW would have been 
smooth, but diverging at the sites of the strongly scattering atoms, each radial part going 
as jv ( Kr R') — tanatRiL' (k) ny (nr R >) . In order to get more feeling for SSWs, let us consider 
some limiting cases: 

If we specify qrl = a — > for all channels, we obtain the bare spherical waves. These are 




Figure 1: (above). Boundary condition 
for the 
screened spherical wave, %fj RL (e, r R ) . Only 
strongly-scattering channels are indicated 
and all hard-sphere radii are equal. 

Figure 2: (right). Screened spherical 
wave centered at the origin, ipQ (e, r) , and 
its slopes, S R0 (top). The same for its first 
energy-derivative function, (e, r) (bot- 
tom). 
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known analytically but, except maybe for small molecules, we never use them. Nevertheless, 
with the normalization specified above they areJl 
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These expressions hold only for r 3> a. Unlike screened spherical waves, the bare ones are 
eigenfunctions of angular momentum and are independent of the surroundings. For positive 
energies they have long range. Like all screened spherical waves, the normalization of the 
bare ones is such that they are dimensionless and, unlike the Bessel and Neumann functions, 
they depend little on energy near the hard spheres. 

Another case is when we specify clrl—O'R for all L, and take e=0. Then ipQ L (0, r) is pro- 
portional to the electrostatic potential from a 2'-pole at the origin, surrounded by grounded 
conducting spheres with radii a# centered at the other sites. Since the hard spheres at 
the neighbors break the spherical symmetry around the origin, the SSW has pure angular- 
momentum character merely at its own sphere, and this holds only as long as the own sphere 
coincides with all other spheres concentric with it. Changing the energy will not change the 
SSW much. If we now let clril' be zero for high Z"s, the SSWs will "wobble" into the hard 
spheres. 

If the hard-sphere radii are generated from repulsive potential wells,! the SSWs are the 
"impurity states" for that repulsive MT potential. 

Since the strongly-scattering components of the SSWs are forced to vanish at all sur- 
rounding hard spheres, the range of the SSWs depends on the choice of hard spheres and 
energy: Consider the spectrum el of the wave equation with the homogeneous boundary 
condition that the solutions vanish at all spheres. This spectrum has a continuum starting 
at s", which in the absence of screening is at zero and which rises with increasing hard-sphere 
radii. Now, the SSWs are localized or delocalized depending on whether their energy is below 
or above the bottom of the continuum. Since we choose energy- independent boundary con- 
ditions for the SSWs, their energy dependence merely enters through the wave equation, that 
is through their curvature, and is therefore small when the wavelength exceeds the diameter 
of the largest interstitial in the hard-sphere solid. 
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If all hard spheres centered on the same site would coincide, then the hard spheres would 
have to be smaller than touching because, if two spheres had a point (or a circle) in common, 
then each one of the SSWs centered on the two spheres would be required to be both zero 
and non-zero at that point (or circle). When only a few low-/ channels scatter strongly, 
neighboring hard spheres may intersect. With decreasing hard-sphere interstitial, the SSW 
sets thus in general become more and more localized, until the hard spheres start to intersect. 
Since from there on, the SSWs are forced to change rapidly near the common circles, their 
behavior becomes chaotic as the circles grow. 

We shall generate the screened spherical waves from the bare ones, because those are 
the only ones we know analytically. Hence, we first consider the question of how to expand 
an arbitrary wave-equation solution, \I> (e, r) , which is regular in all space, except possibly 
at the atomic sites, in an SSW set, {i> RL (e,r R )} , with the same energy. If the number of 
atoms is finite, this energy is supposed to be negative. Moreover, since \1/ (e, r) is a solution 
of the wave equation, the SSW set is supposed to have no weakly-scattering channels and, 
a priory, we treat all channels as strong scatterers. Finally, we shall not truncate the SSWs 
inside the hard spheres, but let them continue to the centers. We now expand \I/ (e, r) in 
spherical harmonics on the hard spheres of the SSW set, thus obtaining the coefficients, 
^ mm (z, dRim) ■ Unless all of these vanish, the linear combination converges to \I> (e, r) : 



because, by construction, the linear combination is a solution of the wave equation with the 
proper energy, and this solution matches \I/ (e, r) channel by channel. In order to convince 
oneself that the latter is sufficient, one may start repeating the argument using an SSW set 
with L-independent hard spheres. In that case, ^ (e, r) coincides with the linear combination 
on a closed boundary, because in the case where the system is infinite such a boundary is 
formed by the entity of all hard spheres, and in the case where the system is finite, the 
boundary is formed by the hard spheres plus the infinity, where both \& (e, r) and the linear 
combination vanish since the energy is negative. 

If all the coefficients ^ mm (e, ^Rim) vanish then \l> (e, r) is an eigenfunction of the hard- 
sphere solid. In this case a complete set must include, in addition to the SSWs, the degenerate 
eigenfunctions, or we may choose a different SSW set for the expansion of * (e, r) . 

Changing the hard-sphere radii, but not the sites and the energy, produces another set 
of SSWs which is also complete in the above-mentioned sense. All such sets are therefore 
linearly dependent. A set of hard-sphere radii is said to specify a representation and the 
transformation from the a to the b representation is obtained by substituting ip R / L i (e, r R >) 
for \1/ (e, r) in the above. Hence, the transformation is 



Ajj I 



lim EE E ^ram (e, r fl) * 'mm (e, a Rlm ) = V (e, r) 



i>\>L> (£, r R') = £ VrL (£, r R) ^RL,R>L> ®Rl) , 



(11) 



RL 



where 




/ (s, o,rl) are the RL components at a^L-spheres of the functions ip b R 



R'U 



(e, r R>) ■ 



11 



Expanding now the left and right-hand sides in spherical harmonics on the 6r» ///-spheres we 
obtain: 

5r»r>5l"L' = Yl ^R"L",rl 0> b R » L n) ip b RL R , L , (e, a RL ) . (12) 



RL 



The two matrices ip a (e, b) and iff (e, a) are thus each others inverses. 

SLOPE AND STRUCTURE MATRICES 

We have specified the SSWs by their nodes and shall need their radial derivatives at the 
hard spheres, that is, the dimensionless slope matrix. Its element S R , L , RL (s) is defined as 
a R 'L> times the //-component of the radial derivative at the a R > £/-sphere, with the positive 
direction taken outwards from R', of ip RL (e, r R ) . This is illustrated in Fig. 0. 

In fact knowledge of the hard spheres and the slope matrix makes generation of the SSW 
set a simple matter: The spherical-harmonics expansion around any site, R', of any member, 
ip RL (e, r R ), of the set is given by radial functions and the function for the L' channel is: 

Vr' L ',rl ( £ > r R') = fv ( £ i a R'L',r R >) 5 R i V . RL + g v (e, a R > L >,r R >) S R , L , RL (e) . (13) 



The local expansion converges for r R i smaller than the distance to the nearest site. In (T3f) , 
fi and gi are solutions of the radial wave equation, [d 2 /dr 2 — I (I + 1) /r 2 + e] rfi = 0, with 
the following boundary conditions for r=a : / has value one and slope zero, and g has value 
zero and slope 1/a. 

The SSW-set, \ip a (e)) , may also be expressed globally as a linear combination of some 
known set, ip b (e)^ , as we saw in Eq. (0). The transformation matrix, ip a (e, b) , is then 
given by Eq. (|T3D with r R i substituted by b R i L i. With the use of Eq. fll3"D , the completeness 
relation (0) thus expresses the transformation from S b (e) to S a (e) . 

In order to generate the slope matrix, we transform to the bare set, which is known 
analytically: Using Eq. fll3l) , S a (e) is expressed in terms of ip a (e, 0) , which is com- 
puted as the inverse of ip° (e, a) . The latter follows from the local, spherical-harmonics 
expansion about R' of the Neuman function centered at R R'): ktii (nr R ) Yl (r R ) = 
T,u3v {™r>) Y L > (r R/ ) B R , V , RL (k) , where 

B R 'l',rl («) = X)47rz- ,+l '- J "CW (« |R - R'l) Y?, tm ,_ n (r^R') (14) 

is the KKR structure matrix, which is Hermitian. The summation runs over /" = 
\V - l\ , \V -l\+2, V + I, and r l+|L! " is real because C LL >L» = SY L (r)Y£,(r)Y L/l (r)dr. 
The on-site elements of B (k) vanish. In this way, we obtain the most important result: 

aS a (e) - aD {j (no)} = [B («) + k cot a (k)]" 1 — ^— , (15) 

j{Ka) ' j(ko) 

where a, j (na) , D {j (na)} , and cota(/t) are diagonal matrices with elements 
orl, {Ka RL ) , -D {j; (Kdfli)} =na RL j' (no) /j (m) , and cot a RL («) . The quantity 
a R'L'S R , L i RL (e) , which is a|j, L , times the L'-component of the radial derivative of ?/>^ L (e) 
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at the //-sphere, form the elements of a matrix which is Hermitian. This matrix, we call 
the structure matrix. 

For the channels to be treated as strongly scattering with the set {ip a } , we take a R L («) 
to be the hard-sphere phase shifts (|3|), and for those to be treated as weakly scattering and, 
thus to be downfolded into the SSWs, we take a R L (k) to be the real phase shifts, T)ri (k) . 
The non-scattering channels do not enter the screening calculation (0), since they neither 
scatter the bare, nor the screened set. The strongly- and weakly-scattering channels thus 
contribute to the size of the matrix to be inverted and the strongly-scattering channels are 
the only ones which will eventually enter the equations for solving Schrodinger's equation. 

Instead of expressing QT3| ) and flT5|) in terms of the usual spherical Bessel and Neumann 
functions, one could of course have divided the factors k, 1 and k' 1 " 1 out on the right-hand 
side of fll^D , or used ipf ( £ , r ) instead of rii («r) , etc.. The only difference between the last 
equation of (|2|) and equation fll^D, is that the Hermitian matrix aS a (e) is normalized in such 
a way as to make its energy dependence as small as possible, and in such a way as to give 
S a (e) a geometrical interpretation, namely as the dimensionless slope matrix. Specifically, 

k tana [B a (k) — K cot a (k)] k~ 1 tana = —j (na) a [S a (e) — D {j (no)}] j (na) , 

so that the screened structure matrices B a (k) and aS a (e) differ because functions of energy 
have been subtracted from the diagonal elements, and because the rows and columns have 
been rescaled with such functions. If we form: 

%>v,rl (e) = -2 (w/a R , L ,) 1 ' [S% L ^ RL (e) + (I + 1) <W, M ] (w/a RL ) l+1 , (16) 

then S a (0) is the conventional {k=0) LMTO structure matrix for the screening constants 

a RL (0) = [2 (2/ + 1)]- 1 (a RL /w) 2l+1 , 

and S a (0) is its first energy derivative for some a (0) . For LMTO users who have developed 
a feeling for the sizes of the conventional structure constants and do not care about the 
new interpretation in terms of logarithmic derivatives, it is of course possible to use the new 
method in the conventional "gauge" (|T6"D. In that case, one must substitute the old poten- 
tial functions, P RL (e) , by -2 (w/a RL ) 2l+1 [D {ip m (e, a RL )} + 1 + 1], with D {ip Ri (e, a RL )} 
evaluated as explained in the following section.! 

Whereas the slope matrix specifies the normal gradients on the hard spheres of all func- 
tions in the SSW set, its first energy derivative, S RLR , L , (e) , specifies the normal gradi- 
ents of the first-energy derivative functions, ip RL (e) , as illustrated at the bottom of Fig. 
0. Since the hard spheres are independent of energy, the energy-derivative functions will 
vanish at all hard spheres, including their own. The first energy derivative of the struc- 
ture matrix in addition gives the overlap matrix of the SSW set: {^ R l (s) I^r'l 1 ( £ )) = 
aS RL RL , (e) . This equation follows from the more general one: {ip R i (e) \^> R 'L' ( £ ')) = 

j (e — e') , which may be derived by use of Green's second 



Srl,r>l' i e ) ~ Srl,r'l> i e ') 
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theorem.B Here, the strongly scattering radial components have been truncated inside the 
corresponding hard spheres as illustrated in Fig.s [1| and |2|, while the remaining, regular 
components extend to the centers of the spheres. 

Considered as functions of e, the eigenvalues of the structure matrix aS a (e) have poles 
when e coincides with an energy eigenvalue, ef, of the hard-sphere solid. For practical 
purposes, the a-radii can be chosen in such a way that the energies of interest to us are 
well below the bottom of the hard-sphere continuum, and below any localized state of the 
hard-sphere solid. For such energies, the eigenvalues of aS a (e) are analytical functions of e. 
This latter point is demonstrated in the left-hand side of Fig. |3], which also demonstrates 
that the energy dependence is weak over the ±10 eV region considered. The right-hand 
side shows that, for low energies or close sphere packings, the slope matrix decays by an 
order of magnitude per shell of neighbors. For a monotonically decaying SSW we expect, 
as illustrated in Fig. |2|, a negative slope at its own hard sphere and positive slopes at the 
neighboring spheres. This is also the behavior found in Fig. ||, at least throughout the first 
three shells. 

In conclusion, the slope matrix generated by inversion of the non singular matrix ([15]) 
contains all the information we shall need about the SSW set. 
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Figure 3: The ssa-element, Sr 00 000 (e) , of the slope matrix for the fee structure and R 
in the O'th to 4'th or 5'th shell. Left: a sp( i=0.7w. Here w is the Wigner-Seitz radius, which 
in the fee structure is 10% larger than that of touching spheres. ew 2 =6.05 corresponds 
to the lowest free-electron energy at the X point. Right: e=0 and a spd =a. The 5-scale 
is logarithmic and the channels with I > 2 were taken as non-scattering. The number of 
atoms in the O'th-5'th shell are respectively 1, 12, 6, 24, 12, and 24. The calculation was 
performed by matrix inversion, Eq. (|15|), in real space for a 79-site cluster. For positive 
energies it was necessary to prevent resonances at the surface of this cluster by enclosing 
it in a concave sphere simulating the boundary condition ?/>=0 on the spheres outside the 
cluster and carrying max(/)=8. The artefact at ew 2 ~6.3 is a surviving resonance. The 
results are accurate only when the SSW is well localized within the cluster. For cases where 
the SSWs are not localized within a cluster of affordable size, we must assume crystalline 
boundary conditions and Bloch-sum B (k) with the Ewald technique before performing the 
screening inversion ([15|) for each k-point. With largely overlapping MT spheres, the energies 
of occupied states are always negative. 



14 



SOLVING SCHRODINGER'S EQUATION WITH KINKED PARTIAL WAVES 

We now come to consider Schrodinger's equation, [—A + V (r) — e,j\ ^ (r) = 0, for a MT 
potential and begin by showing that with our screened spherical waves it is a rather simple 
matter to formulate the matching problem for the solutions, ^ (r) , algebraically: 

First we integrate the radial Schrodinger equation for each strongly scattering channel 
outwards from the origin to the MT radius s R in the potential well v R (tr) , and then inwards 
in zero potential (the MT zero) from sr to the hard-sphere radius orl. The outwards 
integration yields the radial partial wave 4>ri (e, r R ) , and the subsequent inwards integration 
yields the radial partial wave "as seen from free space" ipm (s,tr) , with radial logarithmic 
derivative D {<pm (e, clrl)} at the hard sphere. These two waves match continuously and 
differentiably at sr and they may be seen in the left-hand side of Fig. [|, after multiplication 
by Yl {yr) . Let us assume that <fi RL (e, r R ) and ip RL (e, r R ) have been normalized in such a 
way that <p RL (e, orl) = 1 at the hard sphere; this is what the superscript a here indicates. 
In the case where the hard spheres have been chosen to depend on m, radial functions 
of the same Rl may have different normalizations, hence the subscript L rather than I. 
With this normalization, the free partial wave matches continuously, but with the kink 
Srl,rl ( £ ) ~ D{ip R i (e,a RL )}, to the .RL-projection, ipR IUjRL (e,r R/ ) , of the corresponding 
SSW, ipRL ( £ > t r) ■ Let us furthermore truncate (f) RL (e, r R ) and <p RL (e, r R ) outside the MT 
sphere (0|sij) and, like the SSW, let us truncate ip RL (e,r R ) also inside the a^L-sphere. The 
function [4> RL (s, rR) — ip a RL (e, r R )] Yl (r R ) thus equals the proper partial wave inside the hard 
sphere, where it jumps by — Y L (r R ) , and it vanishes quadratically at the MT sphere with 
a prefactor proportional to the MT discontinuity v R (s R ) . To this function we now add the 
corresponding SSW thus obtaining the kinked partial wave (KPW): 

$ a RL (s, r R ) = [<p RL (e, r R ) - Vrl (e, r R )\ Y L (r R ) + ^rl (e, m) , (17) 

which is also shown in Fig. ^. This function is everywhere continuous, but has kinks of size 
s r'l',rl ( £ ) ~ D {Vri {£, Url)} $r'L',rl at the hard a^-spheres. 

At such a sphere, the kink of the linear combination of KPWs, J2rl ®rl ( e ' t r) c< rl ( £ ) > i s 
therefore Y,rl \Sr'u,rl i e ) - D {Vri Url)} $r>l',rl\ c a RL (e) . If we can now find an energy, 
Si, and coefficients, c RLi , such that 



53 [ S R'L',RL ( £ i) - D ( £ i, a RL}} 5r>L>,RL 

RL 



f RLi = for all R'L', 



then the corresponding linear combination is smooth and therefore solves Schrodinger's equa- 
tion with Si as an energy eigenvalue. 

The statement that the "kink- cancellation condition" ([18j) leads to a solution of 
Schrodinger's equation is exact only for a non-overlapping MT potential. Before continuing 
to the case of overlapping potentials, let us scrutinize our proof a little closer. Each KPW is 
constructed to be a solution of Schrodinger's equation at energy e, except in all shells between 
concentric MT- and hard spheres, and except for the kinks at the hard spheres. In the case 
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where we choose all concentric hard spheres to coincide with the MT sphere (clrl=sr) , 
all shells vanish so the statement is obviously true. That it holds also when the hard 
spheres are different from the concentric MT sphere, follows from the fact that for a linear 
combination with all kinks cancelled, each (p R , L , (e*, r R i) c R i L , i matches the i?' //-projection, 
T.rl ^wl',rl { £ h tr>) c a RLA , of the linear combination of SSWs, J2rl ^rl { £ h r fi) c a RLti , at a R , L , 
in value and in slope, and since both radial functions are solutions of the same second-order 
differential equation, namely the /"th radial wave equation, they must be identical. As a 
consequence, 

<P%i/ ( £ ii r R') c R'v,i ~ ^r'l',rl ( £ h r R >) c a RLi = 0, for < r R > < s R , and R'L' G strong scat. 

RL 

(19) 

Inside the s^'-sphere then, only terms which satisfy Schrodinger's equation remain, 
namely Yy (yr>)J2rl 4>r'l',rl ( £ ii r R') c RL,i an d ^ ne weakly- and non-scattering channels of 
T.RL V&l i £ i,r R ) 

The KPW is defined in ( |TT| ) as the SSW plus the central, pure angular-momentum con- 
tribution (f) — ip, which vanishes quadratically at the MT sphere. In analogy with Slater's 
augmented plane wave (APW), the KPW might have been named an augmented screened 
spherical wave. This analogy is only complete though when all hard spheres coincide with 
their concentric MT sphere. 

Next we consider the case of MT overlap. Suppose that we have solved the kink- 
cancellation equations fll8|) with logarithmic derivatives calculated for potential wells which 
overlap. To what extent is the resulting smooth function, (r) = J2rl &rl ?r) CRL,i, 
a solution of Schrodinger's equation for the superposition of these overlapping wells? The 
situation is sketched in Fig. |5] and the answer is, that the smooth superposition of KPWs 
solves Schrodinger's equation to leading (first) order in the potential overlap. 
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Figure 4: Kinked partial wave (KPW), |$) = \<f>) - \<p) + \ip) , and LMTO, \x) = |$) - 
$) K~ X K, for Si Px+y+z plotted along the [11 Indirection towards a nearest neighbor in the 
diamond structure. Note the change of length scale between the left and right panels. 
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Since we have only considered the strongly-scattering channels in this one-dimensional 
figure, let us now be a bit more careful. Using the definition ( |TTD and the following defi- 
nitions: ipi (r) = Y.rl iPrl ( £ h tr) crla, <P R (e*, r R ) = El <Prl ( e u *r) crl,i, and similarly for 
ip R (e h r R ) , we obtain: 



-A + J2 V R ( r R) ~ £ i 
R 



*i (r) 



E 

R< 



VR>(r R >) E [0a ( £ ii r R) ~ Pr. ( £ ii r n) 

RytR' 



+ E h A + V R M - £ i] Wk ( £ h m) ~ <Pr ( £ i, r R )} + 



R 



R 



fa W 



E V R' i r R') E t^fl - Vr (ei, r R )) -J2 V R ( r R) [<Pr ( £ h v r) ~ fa ( r )] ~ [ A + £ i\ fa ( r ) 

R' R^R' R 



E V R' ( r «') E Wk ( £ " T R) ~ Vr (ei, r fl 

1 



9 ( r «') E vr(s r ) 



j pairs 

~ E ( r «' 

Z RR' 



sr - r^) 2 + o ((s R - r R 

Vr (sr) *j (r 



[sr> - r Rl f + (s R - r R f 



(20) 
(21) 

(22) 





Figure 5: Middle: Overlapping 
potential wells, va (ta) an d vr (tb) , 
centered at sites A (far left) and 
B (far right). Top: The solu- 
tion (f) A = Y,L ( l ) AL(£i,r A )cAL,i joins 
smoothly onto the free solution ip& at 
sa- <Pa runs backwards to a a, where 
the kink with the interstitial solu- 
tion, E^ = Ei?L IpRL ( £ i> r R} c RL,ii 

is cancelled. Similarly for <f)R 
and <ps- Due to kink cancellation, 
the resulting wave function, (ft a — 
<Pa+ 4>b - <Pb + E^, equals <f) A + 
4>b — E^ in this picture where the 
angular-momentum character has been 
suppressed. Bottom: The error, 

[-A + V A + V B - Si] \<j) A + <j) B ~ E fa 

= v B |0a) + v A I 0b) - (va + v A ) \T,fa 
= v B \4>A-J2fa + V A \(pB -J2fa , 
consists of two terms each of which, 
e.g. the first, is the product of 
vb(j"b) and 4>a — EVs which vanishes 
like v A (sa) (sa ~ ^a) 2 4>a (sa) at the 
sa boundary. Hence the error is of sec- 
ond order in the potential overlap. 
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Here, we have first of all made use of the fact that (r) is smooth so that we can apply 
the A-operator to its individual, kinked or discontinuous parts without keeping track of all 
the resulting diverging terms, because they will cancel in the end. In obtaining the 3rd line, 
we have used that <ft solves Schrodinger's equation for its own well. Eq. ( PD| ) has then been 
obtained by use of Eq. flUf ) for the strongly-scattering partial waves, plus the fact that 
the weakly- and non-scattering channels (A) of ipi (r) solve Schrodinger's equation, i. e. that 

[ A + e] 4>%l ( £ > r n) = vr> (r R/ ) £ A 4>%a,rl 0, r i?) • 

Returning to the strongly-scattering channels of J2r v r ( r R,) VPr { e h r i?) — ipi ( r )] > if the 
overlap is so large that the Sij-sphere overlaps a neighboring a r'l> -sphere, then it is simplest 
to imagine that we have not truncated the SSWs inside their hard spheres, because otherwise 
the cancellation (O) would not take place inside the sr-oril' overlap. For consistency then, 
we should not truncate the free partial waves tp inside their own hard sphere either. The 
resulting divergencies at the sites, of the SSWs and of the free partial waves, of course cancel 
for the smooth linear combinations. This undoing of the truncation inside the hard spheres 
is not necessary, but it simplifies the bookkeeping. 

The result (|20p is then in agreement with what we found in Fig. |5|, that the error is a 
function which vanishes outside the regions of overlap and that inside such a region, it is the 
product of a function, vr> (tri) , which vanishes with a small discontinuity at one of the MT 
spheres and a function, 4>r — ^>r, which vanishes quadratically at the surface of the other MT 
sphere, with a prefactor proportional to the discontinuity of that MT potential. Remember 
that the radial part of ifiR is supposed to continue to the origin. The result, which is given in 
(pl|), may be obtained from the radial Schrodinger and wave equations. Finally, in expression 
(B^) we have kept only the term of leading order and have used that <p R tr) ~ ^ (r) in 
the region picked out by the other factors. Hence, the error of the wave function is of second 
order in the potential overlap. 

The error of the one-electron energy may be obtained by first order perturbation theory 
as: A^j = Si — ef uc ~ — |— A + J2r v r ( r i?) ~ £ i \ ^i) an d, to leading order, we find from 
Eq. (p2|) that the error of the band-structure energy is@ 

occ pairs /-p _i_ T»'\ 

~ -^T,^-^\ 5 ^RR'MsR)v RI (sR l )pl [ — r - j, (23) 

where lorri = _^ — 1 is the radial overlap. (24) 

|R. — R- 1 

In a last section we shall demonstrate how this works for the third-generation LMTO method. 
Finally, it may be noted that the appearance of the KPW in Fig. f| is hardly influenced by 
the MT overlap. This figure in fact applies to an overlap of u—14%. 

Like the slope matrix, the kink matrix is not Hermitian, but the matrix 

k r'L',rl (e) = a R'L> [S% L ^R L (e) - D {ipRi (e, a RL )} 5r, v> rl] (25) 

is.0 This matrix is the renormalized screened KKR matrix. If we multiply each of the kink- 
cancellation equations fll9D with the corresponding hard sphere radius, clr/l' j these equations 
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take the form: J2 Kc = and, hence, they are the screened KKR equations. Just as the 
first energy derivative of the structure matrix is the overlap matrix for the set of SSWs, so 
the first energy derivative of the KKR matrix is the overlap matrix for the set of KPWs.0 
In fact, one may show that the KKR matrix itself is the energy minus the MT Hamiltonian 
in the basis of the KPWs with the same energy, that is, 

K*u,rl (e) = (*) k " (- A + V) | <S> RL (e)> . (26) 

For Green-function and CPA calculations it has been very important that the transfor- 
mation (0) of the resolvent, [P (z) — S] _1 , from one representation to another is merely a 
scaling rather than a matrix operation. This turns out to hold also in the new formalism, 
and it means that such calculations may now be performed with more realistic potentials 
and including downfolding. The result is: 

K h (z)" 1 = a~ 1 g a (z, b) <p a {z, b) + <p a {z, b) K a (z)' 1 <p a {z, b) (27) 



and has been obtained by use of the completeness relation (12), the one-centre expansion 
([□D, and the following Wronskian relationsJl 

ag b (a) = -bg a (b) , af b (a) = b 2 g a (b)' , a 2 g b (a)' = bf* (b) , a 2 f b (a)' = -b 2 f a (&)' , 

where the common energy argument, z, has been dropped. 

LOW-ENERGY, FEW-ORBITAL, TB HAMILTONIANS; HTSCs 

If the energy dependence of the renormalized screened KKR matrix is linearized around 
some chosen energy e u , 

K a (e) w K a + (e - e v ) K a = - ($ a |-A + V\ $ a ) + e ($ a | $ a ) , (28) 



then the KKR equations (jTSQ have the form of an algebraic eigenvalue problem. In (p8|) 
and in the following, omission of an energy argument e means that the function is evalu- 
ated at e u . The basis set which, by use of the Raleigh-Ritz variational principle for the MT 
Hamiltonian, gives rise to this problem turns out to be the KPW-set at the fixed energy, 
e v . This follows from Eq. ( ^6|) and is expressed in the second part of Eq. (p8|) . Since 
the off-diagonal elements of the overlap matrix, K, only influence the energy eigenvalues 
to order (£j — e u ) 2 , we may even neglect the non-orthogonality of the KPWs and, for a 
crystal, obtain the correct Fermi surface, ei(k)=eF = £ v , an d the correct group veloci- 
ties, dsi (k) /d\s\ e , by diagonalization of a first-order Hamiltonian whose matrix elements 



are simply: — K^ ^y = —Kf^ ^^ j \j K^ ^K^, L , R , L , . This Hamiltonian is completely 

analogous to h 3 in the ASA, but — K a implicitly contains the integrals over the interstitial 
region and the downfolded channels, and it works to leading order in the overlap of the 
potential wells. The range of —K RLjR i L i in i?-space, and the size of the energy window inside 
which the linear approximation holds, depends on the screening. Crudely speaking, the more 
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strongly-scattering channels included, and the larger their hard spheres chosen without being 
touching, the shorter is the range of the hopping, and the wider is the energy window. 

— K a can be used as the low-energy, few-orbital, single-particle part of correlated 
Hubbard-type Hamiltonians, as we shall now demonstrate for a generic high-temperature 
superconductor (HTSC). We have in the past^ been able to derive such a Hamiltonian for 
YBa2Cu3C>7 using the second-generation LMTO package.^ That procedure, however, re- 
quired a lot of hand- work and much insight, and has proved cumbersome to use in general. 
The new procedure is far more automatic and accurate,^ and has already proved successful 
for the ladder compounds.^ 

The basic structural element of all HTSCs is a CuC>2 layer, which is a quadratic lattice 
with copper at the corners and oxygen halfway between all copper nearest neighbors. In the 
left-hand side of Fig. |6], the copper sites are those which carry either a d x 2_ y 2 or an s orbital, 
and the oxygen sites are those which carry a p x or a p y orbital. Different HTSC materials 
have different stackings of the Cu02 layers with various "insulating" and/or "doping" layers 
between them. Nevertheless, the calculated LDA band structures near what is believed to be 
the Fermi level of optimally doped HTSCs are very similar, and similar to that calculated for 
the simplest possible such material; dimpled CaCu02- In this compound, the Cu02 layers 
are stacked in the z-direction and are separated by calcium, which sits in the hollow between 
the eight coppers of the two neighboring layers. The oxygens in the Cu rows running in the 
x- (y-) direction are dimpled out of the plane by + (— ) 7 degrees. The right-hand side of 
Fig. |6| shows a central Cu02 layer seen from the side, with a d x 2_ y 2 and an s orbital on the 
copper sites and a p z orbital on the oxygen site. On the Cu02 layer above is shown a Cu 
s orbital and on the Cu02 layer below, an O p z orbital. Dimpled CaCu02 is a calculated 
structure,!! a theorists dream which hardly exists in this simple form in nature. Its LDA 
energy bands, which we shall now consider, are nevertheless very similar to those calculated!! 
for YBa2Cu307, one of the only known stoichiometric optimally doped HTSCs. 

At the Fermi level there is only one band per Cu02 layer, and this is the anti-bonding pda 
band formed from the Op x - Cu4 2 -!/ 2 - Op y orbitals. This band is at the top of the 10 eV 
broad Op - Cu d complex consisting of 16 bands, the upper (anti- and non-bonding) part of 
which may be seen in Fig. (a). According to the LDA and the so-called Van Hove scenario 
of HTSC, the Fermi level (zero in the figure) for the optimally doped compounds is very 
close to the saddle-point of the conduction band at (ak x , ak y ) = (tt, 0) . Hybridization with 
the Cu s band, which is 5 eV above, has pushed this saddle-point of the anti-bonding pda 
band down in energy, to a point where it just "straddles off" the top of the anti-bonding pdir 
bands. This makes the structure susceptible to out-of-row movements of oxygen, because 
this will mix a and tt bands. In particular the stable structures of CaCu02 and YL^CuaOy 
have oxygen dimpled seven degrees out of the layer, and this mixes Op z character into the 
conduction band in such a way that its saddle-point at (tt, 0) becomes "extended" that is, 
the dispersion towards (0, 0) becomes proportional to k A , i.e. flat, while in the perpendicular 
direction, towards (tt, tt) , it remains k 2 . The mixing pushes the corresponding pdn band 
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down in energy by about half an eV and leaves the top of other pdrt bands about an eV 
below the Fermi level. We thus see that the orbital character of the conduction band, which 
is the only one we wish to describe, is quite mixed. 

The converged LDA bands are showed in panel (a) of Fig. [7]. For comparison, panels 
(b)-(d) show the bands calculated with various "minimal" LMTO sets, specifically, with 
only the six O p orbitals (b), with only the Cu d x 2_ y 2 orbital (c), and with the six O p 
orbitals plus the Cu s and d x 2_ y 2 orbitals (d). These four calculations all employ the full 
3rd-generation LMTO formalism, to be described in the following section, in which the 
Hamiltonian and overlap matrices, (|3T| ) and (0), are given in terms of K a (ef) and its first 
three energy derivatives. Panel (b) and (c) demonstrate the power of downfolding in the 
3rd-generation LMTO scheme: One may for instance completely leave out the Cu d x 2_ y 2 
LMTOs by attaching that partial-wave character to the tails of the neighboring O p LMTOs 
(b), or one may completely leave out the O p LMTOs, keeping per cell just the one Cu d x 2_ y 2 
LMTO whose tail then incorporates the O p, Cu s, and other characters (c). As one can 
imagine, such massive downfolding leads to long range of the LMTOs. As an example, the 
Fourier transform of the conduction band shown in panel (c) is the two-centre Hamiltonian 
in the representation of orthogonalized Cu d x 2_ y 2 LMTOs, where the cone-like feature of the 
band around (0,0) , caused by near degeneracy of the C\xd x 2_ y 2 and Op x orbital energies, 
gives rise to very long range. This long-ranged, single-band Hamiltonian, we have called (the 
single-particle part of) the "physical" low-energy Hamiltonian.il 

What we shall be interested in here is a "chemical" Hamiltonian, which has short range 
and whose TB parameters behave in a meaningful way when the structure is deformed and 




Figure 6: The eight orbitals [O x p x , y ,z, O y p x , y ,z,, Cu s, and Cu d x 2_ y 2] and the values 
of their energies (with respect to E x 2_ y 2) and two-centre hopping integrals (eV). These 
values were obtained as the matrix elements of — K a (ef) with all channels other than the 
eight downfolded. 
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when we proceed to similar materials. Which orbitals such a chemical Hamiltonian should 
contain is then dictated by the range of the corresponding K a (e) matrix. If we imagine 
a Taylor series like (p8"D, it is conceivable that the higher energy- derivative matrices have 
longer range. We therefore expect to obtain the shortest range when the energy region of 
interest is so small that we only need Kr L B!L , (ep) as defined above. For dimpled CaCu02, 
the chemical basis set turns out to be the one used to generate the bands shown in panel (d). 
For the same eight orbitals, we show in panel (e) the bands calculated by diagonalization of 
the effective two-center Hamiltonian —K (ef) ■ We see that this approximation conserves the 
shape of the conduction band in the relevant range of energy. All computations illustrated 




Figure 7: LDA energy bands for dimpled CaCu02 calculated with the 3rd-generation 
LMTO method using a converged LMTO basis, (a), and five different simplifications, 
(b)-(f). Reciprocal-space distances are in units of the reciprocal of the Cu-Cu dis- 
tance. ef= =e u . (a): All but the Cus and d, and the Op orbitals downfolded. 
acus=acud=0.87t, ao p =0.75t, where the touching-sphere radius, t, is 1/4 the Cu-Cu 
distance, (b): All but the six oxygen p orbitals downfolded. ao p =0.96t. (c): All but the 
single Cud x 2_ y 2 orbital downfolded. a Cux 2_ y 2=0.62t. (d): All but the six oxygen orbitals, 
the Cu s, and the Cu<i x 2_ y 2 orbitals downfolded. a Cus =0.87t, a Cu:r 2_ y 2=0.68t, a Op =0.96t, 
(e): Like (d), but with H — ef = —K and O = 1. (f): Like (e), but with K truncated after 
3rd-nearest-neighbor hoppings, that is, with the orbital energies and two-centre hopping 
integrals given in Fig. [| 
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so far were converged in R-space. In fact, they were performed in k-space, which means 
that we started out using the Ewald method to compute B (re, k) . When we now Fourier 
transform Kf lLR , L/ (ep,k) , we find that the only non- negligible matrix elements are those 
given by the orbital energies and two-center hopping integrals in Fig. []. Panel (f ) of Fig. |7| 
shows the corresponding TB energy bands. This orthogonal, two-center TB Hamiltonian, is 
seen to reproduce the conduction band very well and to give a satisfactory description of the 
neighboring bands. This TB Hamiltonian, which we have generated almost automatically, 
could also have been calculated without the Ewald scheme, by inversion of Eq. (|T5| ) in It- 
space. A bit of trial and error is still needed in finding an optimal choice of the hard-sphere 
radii. The ones we used are listed in the figure caption. 

LINEAR MUFFIN-TIN ORBITALS 

The first-order Hamiltonian —K does not suffice to describe the energy spectrum over 
the 10-20 eV range spanned by the valence and lower conduction bands of strongly bonded 
materials. Nor does inclusion of terms beyond the linear in the Taylor series help, 
because this does not lead to an algebraic eigenvalue problem. What is needed, is a set 
of energy independent orbitals which, in contrast to the set of KPWs at a fixed energy, is 
complete to linear order in e — e v . 

From a set of KPWs, we first define a set of energy dependent MTOs: 

\ X (e))= {^(e^-^k-'Kie) (29) 

Here and in the following we often drop the common superscript a, and omission of an energy 
argument means that e=e v . Moreover, we have used the notation in which \x{ £ )) is a row 
vector with elements \xrl (e)) = Xrl (e, i\r) and K is a matrix. <&rl (fr) is the first energy 
derivative at e v of the KPW, §rl (e, Fr) , defined in fllTI) . Since the hard spheres are kept 
independent of energy, the strongly-scattering channels of the energy- derivative functions $ 
vanish at all the hard spheres. The ^-part is sketched in the bottom half of Fig. |2] and the 
Si p x +y+ z MTO at energy e v , that is the LMTO, is shown together with the corresponding 
KPW in the right-hand side of Fig. |j. 

The superposition of ^-functions added to the KPW in (P9|) is such as to make the 
MTO smooth. That this is so is seen immediately by forming the kink matrix for the MTO: 
K (e) — KK~ l K (e) = 0. Still, the set of MTOs remains complete with respect to the MT 
potential, because with being the energy and Cj a corresponding solution of the KKR 
equations, K (e t -) c» = 0, we find that the same linear combination of MTOs is: \x (£*)) Cj = 
|<3> (gj)) Cj = l^j) . In contrast to the KPW, the MTO is independent of energy to linear order 
because by differentiation of ( |29|) with respect to energy and subsequent setting e=e v we get: 
\ x ) = $^ - <j>\ K~ 1 K = 0. The energy-independent set of LMTOs, |x) = |$) - $) K^K, 
is therefore complete to linear order with respect to the MT Hamiltonian and therefore 
yields eigenvalues with errors proportional to (e$ — e u ) A . For comparison the conventional 
single-K LMTO set is complete to zeroth order in the MT interstitial, albeit to first order in 
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the spheres, and therefore yields eigenvalue errors of order (si — e„) which originate from 
the interstitial. This is illustrated in Fig. 8. A price for carrying not only ip, but also ip 
functions, is that the new LMTO sets corresponding to different hard-sphere radii are no 
longer linear combinations of each other; the wave-function error, A a ■ (ei — e u ) 2 , has an 
a-dependent pref actor. 

We now derive the expressions for the Hamiltonian and overlap matrices in the new 
LMTO basis. For the integrals in all space of KPWs and their first energy derivative func- 

K • The 



tions, one obtains: ($ | <3>) = K, | $ 
LMTO overlap matrix is therefore: 



$ I $ 



^K, and 



$ I $ 



3! 



(x I x) 



($ I $) - ($ I $\ R- L K - KR- 1 ($ | $) + JCfiT 1 | $) K~ X K 
= if - 1 [kk- x K + KK~ X K) + ^i^" 1 # if" 1 ^. (30) 

The matrix elements of the MT Hamiltonian used to generate the LMTO set may be found 
in a similar way. Since the LMTO is smooth there are no problems with Hermiticity like 
those occurring for the matrix elements between KPWs alone. What we mean is, that 
the result (^6|) cannot be obtained by naively taking matrix elements of an equation like: 



[H — e] |$ (e)) = 0, where H = —A + V, or of its energy derivative: [H — e u \ $y = |$) . For 
matrix elements between smooth linear combinations of KPWs like: 

{x\-A + V -e v \x) = ($|#-e„|$) - (§\H - e u \§) K~ l K - KK~ l (§\H - e u \<$> 
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Figure 8: LMTO errors of the lowest free-electron energies (0, 0.75, and 1 Ry) at the 
fee T, L, and X points as functions of the energy-expansion parameter e v . The states 
at L and X are doubly degenerate and split when e v ^ e. (a): Old LMTO method with 
sprf-basis and s—w. (b): New LMTO method, Eq.s ( ffiJj) and fl3"T]) , with spd-b&sis and 
a sp d=0.7w=0.77t; the s-value is irrelevant. In (a) and (b) the bare structure matrix was 
Bloch-summed with the Ewald technique before it was screened by matrix inversion, Eq. 
(jjD, in k-space. (c): Like (b) except that the inversion dl5|) was performed in R-space 
using a 79-site cluster enclosed in a concave sphere. 
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+KK~ l ($ \H - e v \ $) K^K 

-($\H-e v \ $) i^if + KK- 1 ($ |# - e„| 6) 

- ($ | $) K^K + KK^ 1 ($ | $) K^K 



= -K + -KK- 1 kk- 1 K (31) 

such procedures are however correct when used consistently for all terms. Expression ( |3~TD 
thus gives the MT Hamiltonian matrix which, together with the overlap matrix (^), are 
given exclusively in terms of K, K, K, and K ■ These matrices are square and labelled by 
the channels of the strong scatterers. We stress, that in the 3rd-generation LMTO method, 
downfolding takes place at the screening stage (TT5|), where it removes the weakly-scattering 
channels from the structure matrix S (e) . The calculations for CaCu02 presented in Fig. |7| 
(a)-(d) employed this formalism and convincingly demonstrated the new downfolding. 

An approximation, which goes beyond the ASA and is not based on dividing space into 
spheres and neglecting the remainder, consists of neglecting all off-diagonal elements in the 
rea/-space representation of K, K, and K. With this new ASA, we have avoided the matrix 
inversion, K~ x , and the formalism contains only one matrix, which we may take to be the 
first-order two-centre Hamiltonian —K defined in the previous section. This corresponds to 



renormalizing each KPW and each MTO according to: (e)> = \Qrl (e)) / y (®rl) 



®rl (e)) / <JK RL)RL and \\rl (e)) = \xrl (e)) / V K RL)RL , and the rows and columns of 



the KKR matrix accordingly: K RL)KV (e) = K RL ^ R , V (e) j y 1 K RL . RL K R , L ,, R , L > . With this 

renormalization, and taking e=e u , it is easy to see that expressions ( |30l) and (|3l|) reduce to 
the simple ASA form (Rp and (^). 

We can develop an exact formalism by Lowdin orthonormalizing the KPWs, instead of 
merely normalizing them: The overlap matrix for the renormalized KPWs is: 

^RL | &R>L>") =k R L, R 'L'= 3rL,R'L' + ArL,R'L', (32) 

where A is a Hermitian matrix with vanishing diagonal in i?L-space. Its off-site elements 
(R 7^ R') are usually considerably smaller than unity and if we now define a Hermitian 

i.-l/2 ' _i 

matrix: K =(1 + A) 2 =1 — |A + |A 2 — which is the power-series expansion in 

- \ ~ \ ~~ 1/2 

A, then the linear combinations = $(e)> K are seen to form an orthonormal 

set when e=e u . This is formally like in the conventional ASA. The partial waves truncated 
outside and normalized inside the atomic s-spheres become in the formalism of the 3rd 
generation the Lowdin orthonormalized kinked partial waves. The transformed MTO set is: 



■ -1/2 



|x(e)>= \x RL (e))k = |<i(e)) + 



$ ) h (e) , where 



~ - 1 / 2 - ~ - 1 / 2 / I 3 \ ~ / 1 3 

/, (.:) = ~K K(e)K = ( 1 - -A + -A 2 - ...\K{e) (l - - A + -A 2 - 



(33) 
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Since this expression for the MTO set is also formally identical with an expression which, 
with the old definitions, was valid only in the ASA, everything else works out the same. E.g., 



we find: 



X 



$> + 



$^ h = 0, because h = — 1. The Hamiltonian and overlap matrices 



are thus given by (g) with h = h (e u ) 



$|$^ = ($|<^ = -_ an dp + ^ = -^|$^ = -^|$^ = ^|$^ = -ll. (34) 

In the 3rd-generation LMTO, h, o, and p are square matrices labelled by the strongly- 
scattering channels. What we have accomplished is therefore to transform the new Hamilto- 
nian and overlap matrices, (|31|) and (0), into the form (|8|), which was previously valid only 
in the ASA. 

In this language the new ASA corresponds to neglecting A as well as the off-diagonal, 
real-space parts of o and p. A better approximation is to keep A to first order in Eq. (|33|) , 
and then to neglect the off-diagonal parts in the real-space representation of o and p. In 
this way we still need to specify only one matrix, namely the first-order, two-center TB 
Hamiltonian, h, at the expense of increasing its real-space range somewhat beyond that of 
—K. In the full formalism we have to specify 2 matrices, the Hamiltonian and the overlap 
matrix or worse, the 3 matrices: h, ( h — — l) , h\ and h, or even worse, the 4 matrices: 
K, K, K, and K whose real-space range increases with the number of energy derivatives 
taken, that is, in order of decreasing importance for the bands near e v . 

Some of this is illustrated in Fig. |9] where we compare the LDA band structure obtained 
from a converged 3rd-generation LMTO calculation (full line) with results (dashed lines) 
obtained using various minimal basis sets, sp 3 in (a)-(c) and sp 3 d 5 in (d), and various trun- 
cations. The empty-sphere spd- and, in (a)-(c), the Si ci-channels were downfolded. Here 
panel (a) demonstrates that it is possible with merely an sp 3 set to obtain an accurate first- 
principles description of the valence and four lowest conduction bands, provided that we 
allow the set to be so long ranged that its Hamiltonian and overlap matrices, (|3l|) and (|30|) , 
extend to 12th-nearest neighbors. This basis is defined by: a s =l.lt, a p =1.0t, and e u = — 2 
eV. As usual, t is half the nearest-neighbor distance. If an accurate sp 3 TB-description is 
needed of merely the valence band, then it is possible to limit the range of the orbitals to 
the extent that the Hamiltonian and overlap matrices can be truncated after the 6th-nearest 
neighbors. In (b) this is achieved mainly by shifting e v down to the middle of the valence 
band. In (c) and (d) we have simplified the calculation of the Hamiltonian and overlap 
matrices by evaluating ([33|) to only first order in A, and by neglecting the off-diagonal ele- 
ments in i?-space of o and p. As mentioned above, this also makes it necessary to tabulate 
only one two-centre matrix, h. (Note that the screened two-centre matrices cannot be com- 
pletely specified by Slater-Koster two-centre integrals like (||), because the screened KPWs 
and LMTOs do not have pure angular-momentum character). Comparison of the dashed 
lines in (b) and (c) shows that this simplification works for the valence-band structure, but 
that the quality of the conduction band, which was not aimed at here, has deteriorated. 
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So far we have not been able with our first-principles procedure to find parameters which 
will decrease the range of the sp 3 first-order two-centre Hamiltonian, h, below 6th-nearest 
neighbors. However, with an sp 3 d 5 basis this is possible, because then also the (i-channels 
can be used for screening. This is demonstrated in panel (d), where the sp 3 d 5 -set with the 
parameters a s =a p =1.0t, aa=0.9t, and e u = — 6 eV, plus the above-mentioned simplification, 
yields an h which can be truncated after 3rd-nearest neighbors. The resulting valence band 
is good and the conduction band very reasonable. 

In the past there have been several attempts to model the energy bands of Si by a simple 
TB Hamiltonian and the need for TB total-energy representations to provide inter-atomic 
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Figure 9: LDA energy bands of diamond structured Si in full lines and various TB approx- 
imations thereto in dashed lines (a)-(d). The corresponding LMTO sets and the real-space 
truncation of the Hamiltonian and overlap matrices ( pT|) and ( pOf ) in (a) and (b), and of h 
(|33|) in (c) and (d), are specified at the bottom of the panels. 



27 



forces for molecular-dynamics simulations has renewed this interest. These attempts range 
from simple nearest-neighbor, orthogonal parametrization of diamond structured Si with an 
sp 3 basis in the 7O's0 to recent work with long-ranged non-orthogonal sp 3 d 5 basis setsS with 
a hope to provide transferable parameters. All these works relied on fittings of energy bands 
and total energies obtained from first-principles calculations. Our method is free from such 
fitting procedures and is purely deterministic. The recent work of McMahan and Klepeisil 
is more similar in spirit to ours, but being based on a full-potential multiple-kappa LMTO 
calculation with the need for subsequent contraction to a minimal sp 3 d 5 basis set, it is more 
complicated and computationally far more demanding. In fact, our method is so fast, that 
for us, transferability is no issue. But in all fairness, our total-energy and force calculation 
is still pending. 

GETTING RID OF THE EMPTY SPHERES 

The full LDA potential for diamond-structured Si is shown in the top left of Fig. |TD|. What 
was used in the LMTO calculations of Fig. however, was the conventional ASA potential 
shown in the top-right panel of Fig. [10], which is slightly overlapping [u;=14%; see Eq. (|24|)1 
and, in addition to the Si-wells, has repulsive wells at the E-sites to describe the hills of 
the potential. Despite its crude appearance, this ASA SiE-potential, gives nearly exact LDA 





Figure 10: Full LDA potential (in Ry) and various MT approximations for diamond- 
structured Si in the (llO)-plane. The two pseudo potentials at the bottom are least- 
squares fits to the ASA potential. 
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valence and conduction bands. But this is a special case. In general, the potential consists of 
spherically symmetric craters with hills in between, and the latter can be of any shape. Such 
a potential is naturally modelled by a superposition of atom-centered spherically-symmetric 
wells, and since we have proved in Eq. (EUp that the KKR method can handle such a 
potential, unless the overlap is too large. The questions are whether this holds also for the 
new LMTO method, and whether the overlap allowed by these methods is sufficiently large 
that the MT zero moves up close to the hill tops and the wave functions tail properly off into 
the voids. Non-MT perturbations would then be local and simple to include. Therefore we 
first try to treat diamond-structured Si. Two appropriate potentials with respectively 30% 
and 60% radial overlap are shown in the bottom panels of Fig. [It]. 

We thus want to fit the full potential, V (r) , to a constant (the MT zero) plus a superpo- 
sition of spherical wells: V (r) ~ V mtz + J2 v r ( r J?) = V (r) . If we decide on a least-squares 
fit, that is, minimization of [V — V] 2 , then variation of the functions Vr (tr) leads to a set 
of coupled integral equations, one for each R saying that the spherical average around site 
R for radius tr should be the same for the sum of the MT wells as for the full potential 
minus V mtz . Variation of V mtz leads to one equation saying that the average of the MT and 
the full potential should be the same. These equations are fairly simple to solve numerically, 
but they do not quite express what we want, because a volume element in a region like a 
void, where the electron has little chance of being, enters with the same weight in the fitting 
as a volume element in say the bond region. What we really want is a pseudo potential 
which, for a certain band, say the valence band, minimizes the mean squared deviation of 
the one-electron energies, Trp [H — H] 2 = Trp [V — V] 2 , and this then brings in the electron 
density, p (r) , as weighting function. This weighting presents little problem for the SV mtz - 
equation, which is merely: / [V — V] pd 3 r = 0, but it complicates the Svr (r^)-equations so 
much, that we decided on keeping p in the ^V^-equation only. Our MT pseudo potentially 
thus pseudizes the hills rather than the core regions. 

Since at this stage, we merely want to see whether we can get rid of the empty spheres 
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Figure 11: Radial behavior of the 
Si and E potential wells of the ASA 
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in the diamond structure by comparing the valence-band structure calculated for the true 
potential with that calculated for its pseudo potential, we take the true potential to be one 
for which we can solve Schrodinger's equation with high accuracy, namely the ASA potential 
shown in the upper right of Fig. [T0| Since this potential is discontinuous at the surfaces of 
the Si and E spheres, its pseudo potentials, shown at the bottom, are not only discontinuous 
at s, but also at the Si AS radius. The radial behaviors of the Si and E wells of the ASA 
potential, as well as that of the Si pseudo potential with 40% radial overlap, are shown in 
Fig. Illl By comparison of the pseudo potentials with 30% and 60% radial overlap shown 



at the bottom of Fig. [10], it is obvious that the latter resembles the true potential most 
closely. Whereas the MT zero of the 14% overlapping ASA potential is only slightly above 
the bottom of the valence band, that of the 40% overlapping pseudo potential lies 6 eV 
higher, and that of the 60% overlapping potential is at the top of the valence band. 

We have now used the new LMTO method [Eq.s fl30|) and (0) with a Si sp 3 d 5 LMTO set 
and the Si /-channels downfolded] to calculate the energy bands for the valence-band pseudo 
potentials as a function of the radial overlap u. The rms and mean errors of the calculated 
valence bands are shown in Fig. [12] by diamonds. Since for increasing u, the potential has 
increasing range and, hence, increasing freedom, the rms error initially falls, but it eventually 
rises again as the kinetic-energy errors given by Eq. (|20| ) and proportional to u> A take over. 
The minimum rms error of 80 meV per electron is reached at 30% overlap. The mean error 
we had expected to vanish for overlaps so small that the kinetic-energy errors are negligible, 
because the pseudo potential was constructed such that / [V — V] pd 3 r = 0. Nevertheless, 
the computation yields a "background" mean error of -50 meV per electron. This is most 
likely due to errors of second order in V — V caused by the unphysical discontinuities at the 
E-spheres of the ASA potential. We expect this background error to vanish and thereby the 




Figure 12: Rms and mean errors of the valence-band energies arising by pseudizing of the 
ASA potential of diamond structured Si in order to get rid of the E wells. The abscissa is 
the radial overlap, u>, as defined in Eq. Q2"4j). For space-filling Si+E spheres, cj=14%, and 
for Si spheres 43%. Overlap correction I modifies the pseudo potential, while II includes the 
proper kinetic energy in the LMTO Hamiltonian. 



30 



rms error to be reduced, when for V we use the full potential in the top left panel of Fig. |10[ 
Since the kinetic-energy error is negative, it represents an attraction between overlap- 
ping atoms, and this might cause problems in molecular-dynamics calculations. However, 
although this attraction increases rapidly with overlap, it does decrease for decreasing inter- 



atomic distance and fixed s-radii [see Eq. (|23|)1. 

If radial overlaps in excess of ~30 % are needed, then the kinetic-energy error must be 
corrected. We have tried two schemes, the results of which are given in Fig. [12] by the stars 
and the triangles. In the first scheme (stars) we have merely modified the pseudo potentials 
by including in the SV m t z -equation the kinetic-energy error to leading order as given by Eq. 
([23]), whereby this equation becomes: / [V — V]pd 3 r = yjj A£j. This leads to a reduction 
of the overlap error, mainly through reduction of the discontinuity v (s) . This correction is 
very simple, but as seen from the figure, hardly sufficient because it only treats the error 
proportional to v (s) 2 uA Our work on the second scheme (triangles) is still in progress.^ 
Here, we evaluate the LMTO Hamiltonian matrix properly to all orders in the overlap, that 
is, we calculate the LMTO matrix elements following Eq. (Ell). Of course, this adds terms 



to expression ([!!]) for the Hamiltonian and spoils the beauty of Eq.s (|8|), (|33"D, and (|34T), but 
we wish to prove that we can control the overlap errors of the new LMTO method, and we 
want to investigate how large overlaps we can handle. The preliminary results shown in Fig. 
12] are encouraging. 
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